The contours from the MCMC seem to be too large. I'm going to take some points on the chain and plot the emulator prediciton along with the "truth" at that point and see if they make sense. Additionally, part of my concern is that the errors for the emulator are not right. If I draw a lot of samples from the emulator at that point vs several repops, they should be similar.
In [2]:
from pearce.emulator import OriginalRecipe, ExtraCrispy
from pearce.mocks import cat_dict
import numpy as np
from os import path
In [3]:
import matplotlib
#matplotlib.use('Agg')
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
from chainconsumer import ChainConsumer
In [4]:
training_dir = '/u/ki/swmclau2/des/PearceLHC_wp_z_test_error'
em_method = 'gp'
split_method = 'random'
In [5]:
fixed_params = {'z':0.0}#, 'r':0.18477483}
In [6]:
emu = OriginalRecipe(training_dir, method = em_method, fixed_params=fixed_params)
In [7]:
emu._ordered_params
Out[7]:
In [27]:
float(path.dirname(path.join(training_dir, 'a_1.00000/')).split('/')[-1][2:])
Out[27]:
In [10]:
from collections import OrderedDict
x = OrderedDict([('logMmin', (11.7, 12.5)),
('sigma_logM', (0.2, 0.7)),
('logM0', (10, 13)),
('logM1', (13.1, 14.3)),
('alpha', (0.75, 1.25)),
('f_c', (0.1, 0.5)),
('r', (0.093735900000000011, 34.082921444999997)),
('z', (0.0, 0.5))])
print x
In [26]:
emulation_point = [('f_c', 0.233), ('logM0', 12.0), ('sigma_logM', 0.533),
('alpha', 1.083),('logM1', 13.5), ('logMmin', 12.233)]
em_params = dict(emulation_point)
em_params.update(fixed_params)
del em_params['z']
param_names = em_params.keys()
In [27]:
rp_bins = list(np.logspace(-1,1.5,19) )
rp_bins.pop(1)
rp_bins = np.array(rp_bins)
rpoints = (rp_bins[1:]+rp_bins[:-1])/2.0
In [28]:
from SloppyJoes import lazy_wrapper
In [29]:
# move these outside? hm.
def nll(p):
# Update the kernel parameters and compute the likelihood.
# params are log(a) and log(m)
emu._emulator.kernel[:] = p
ll = emu._emulator.lnlikelihood(emu.y, quiet=True)
# The scipy optimizer doesn't play well with infinities.
return -ll if np.isfinite(ll) else 1e25
# And the gradient of the objective function.
def grad_nll(p):
# Update the kernel parameters and compute the likelihood.
emu._emulator.kernel[:] = p
return -emu._emulator.grad_lnlikelihood(emu.y, quiet=True)
In [30]:
subsample_idxs = np.random.choice(emu.y.shape[0], size = 100)
In [31]:
gp = emu._emulator
gp.compute(emu.x[subsample_idxs, :], emu.yerr[subsample_idxs])
In [32]:
def resids(p, gp, y):
gp.kernel[:] = p
gp.recompute()
return gp.predict(y, gp._x, mean_only=True)-y
In [33]:
print resids(np.ones_like(gp.kernel.vector), gp, emu.y[subsample_idxs])
In [ ]:
results = []
for i in xrange(10000):
if i%1000==0:
print i
vals = np.random.rand(gp.kernel.vector.shape[0])*2
args = (gp, emu.y[subsample_idxs])
try:
result = lazy_wrapper(resids, vals, func_args = args, print_level = 0,h=0.1,\
artol = 1e-9, xrtol = 1e-21, xtol=1e-20, gtol = 1e-9)
except:
continue
results.append(result)
results = np.array(results)
In [ ]:
for i in xrange(results.shape[1]):
plt.subplot(4,2,i+1)
plt.hist(results[:,i]);
In [ ]:
#from SloppyJoes import fdjac
In [ ]:
def fdjac(x, fvec, func, eps, center_diff):
epsmach = np.finfo(float).eps
dx = np.zeros_like(x)
fjac = []
if center_diff:
for i in xrange(x.shape[0]):#TODO vectorize
h = eps*x[i]
h = eps if h < epsmach else h
dx[i] = 0.5
temp1 = func(x+dx)
temp2 = func(x-dx)
print temp1- temp2
fjac.append((temp1-temp2)/h)
dx[i] = 0.0
else:
for i in xrange(x.shape[0]):
h = eps *abs(x[i])
h = eps if h < epsmach else h
dx[i] = h
print dx
temp1 = func(x+dx)
dx[i] = 0.0
print temp1, fvec
fjac.append( (temp1-fvec)/h)
return np.stack(fjac).T #not sure bout the dimension here
In [ ]:
f = lambda x : resids(x, *args)
fdjac(vals, resids(vals, *args),f, 0.1, False)
In [ ]:
for i in xrange(10):
print resids(i*vals, *args)
In [ ]:
emu._emulator.kernel[:] = result
emu._emulator.recompute()